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Abstract 

We present a numerical framework to treat infinite time horizon spatially distributed optimal 
control problems via the associated canonical system derived by Pontryagin’s Maximum Principle. 

The basic idea is to consider the canonical system in two steps. First we perform a bifurcation 
analysis of canonical steady states using the continuation and bifurcation package pde2path, yield¬ 
ing a number of so called flat and patterned canonical steady states. In a second step we link 
pde2path to the two point boundary value problem solver TOM to study time dependent canonical 
system paths to steady states having the so called saddle point property. As an example we consider 
a shallow lake model with diffusion. 
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1 Introduction 

In [BX08] the authors consider economically motivated deterministic optimal control (OC) problems 
with an infinite time horizon and a continuum of spatial sites over which the state variable can 
diffuse. Using Pontryagin’s Maximum Principle (see §2.2 for general background and references on 
OC problems in a PDE setting) they derive the associated canonical system and show the remarkable 
result that under certain conditions on the Hamiltonian there occurs a Turing like bifurcation from 
flat to patterned steady states of the canonical system, and call this phenomenon optimal diffusion- 
induced instability (ODI). Here we present a numerical framework to (a) study such ODI bifurcations 
of canonical steady states (CSS) numerically in a simple way, and (b) study their optimality by 
calculating and evaluating time-dependent paths to and from such CSS. As an example we use one of 
the three examples presented in [BX08], namely a version of the, in the field of ecological economics, 
well-known shallow lake optimal control (SLOC) model, cf. [Sch98, MXdZ03, CB04]. 

We use the acronyms FCSS and PCSS for flat and patterned canonical steady states, and similarly 
FOSS and POSS for optimal canonical steady states, and summarize FOSS and POSS as OSS. The 
SLOC model has up to three (branches of) FCSS in relevant parameter regimes, and in these regimes 
we also find a large number of (branches of) PCSS. In this situation of multiple CSS a local stability 
analysis at a given CSS is in general not sufficient. Here, local stability analysis means that the 
stationary canonical system is analyzed, analogous to a steady-state analysis of the canonical system 
derived from an optimal control problem without spatial diffusion. It is well known and shown for many 
models that the appearance of multiple CSS (even if these steady-states are saddle-points) does not 
necessarily imply the appearance of multiple steady-states in the optimal system, cf. [GCF+08, KW10]. 
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The reason is that there can exist a non-constant extremal solution, i.e. a canonical path connecting 
the state values of a given CSS to some other CSS, and yielding a higher objective value. Therefore, 
to study whether a CSS corresponds to an OSS, the values of the associated stable paths also have to 
be considered. 

Here we numerically compute the bifurcation behavior of FCSS and PCSS for the SLOC model in 
some detail, and study their optimality by evaluation of their objective values J and comparison to 
time-dependent canonical paths. Such a global analysis is inevitable and has to accompany the local 
stability analysis. Since in general the pertinent ODEs or PDEs cannot be solved analytically we have 
to use numerical methods for the calculation of FCSS and in particular PCSS, and for the calculation 
of t-dependent canonical paths. For the steady state problem we use the continuation and bifurcation 
software pde2path [ VR14], based on a spatial finite element method (FEM) discretization, which 

we then combine with the boundary value problem (BVP) solver TOM to obtain canonical paths. 

A standard reference on ecological economics or “Bioeconomics” is [Cla90], which also contains 
a very readable account, and applications, of Pontryagin’s Maximum Principle in the context of 
ODE models, while [GCF+08] focuses more on socio-economical ODE model applications. Besides in 
[BX08], and in [BX10], PDE models roughly similar to our diffusive SLOC model are considered in, 
e.g., [ XV07, AAC11, DHM12, ACKLT13, ADS14], partly including numerical simulations. However, 
these works are in a finite time horizon setting, and with control constraints, which altogether gives a 
rather different setting from the one considered here. See Remark 2.4 for further comments. 

In §2 we present the SLOC model. To give some background, in §2.1 we briefly present the OD 
(ODE) version, and some basic concepts of optimal control, in particular Pontryagin’s Maximum 
Principle. In §2.2 we turn to the distributed case, explain the associated canonical system, and relate 
our method to other approaches to PDE OC problems. In §3 we explain the numerics to first compute 
the bifurcation diagram of canonical steady states, and then to solve the BVP in t. We mostly focus 
on one spatial dimension (ID), but also give a short outlook on the 2D case. It turns out that in 
the parameter regimes studied here the PCSS are not optimal, but nevertheless they play a relevant 
role. Moreover, calculating optimal canonical paths to FCSS yields interesting and to some extent 
counter-intuitive information about the optimal control of the distributed SLOC model. 

In §4 we close with a short summary and discussion, and Appendix A contains remarks about the 
saddle point property for CSS in a PDE setting, starting on the discretized level. 

Our software, including demo files and a manual to run some of the simulations in this paper, can 
be downloaded from www.staff.uni-oldenburg.de/hannes.uecker/pde2path. In fact, the present 
paper is the first in a series of four related works, the other three being [Uecl5b, Gral5, Uecl5a]: 
[ Jecl5b] contains a Quickstart guide and implementation details of the add on package p2poc to 
pde2path used for the computations in this paper. Thus, the reader interested in these details should 
read (parts of) [Uecl5b] in parallel. Next, [Gral5] explains the usage of OCMat orcos. tuwien. ac. at/ 
research/ocmat_software/ to study ID distributed OC problems based on spatial finite difference 
approximations, with the same SLOC model as in the present paper as an example, and thus obtaining 
comparable results, but also studying a second parameter regime. Finally, in [Uecl5a] we apply p2poc 
to an OC problem for a reaction-diffusion system modeling a vegetation-water-grazing interaction. In 
contrast to the SLOC model studied here, this yields dominant patterned optimal steady states in 
wide parameter regimes, and thus interesting new results on spatial patterns in optimal harvesting. 

2 The model, and background from optimal control 

2.1 The shallow lake model without diffusion 

A well known non-distributed or OD version of the SLOC model, see e.g. [Wag03], can be formulated 
in dimensionless form as 

P oo 

J(P 0 ,k(-))-= e~ rt J c (P(t),k(t))dt, (la) 
Jo 


V(P 0 ) := max J(P 0 , &(•)), 

fc(-) 
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(lb) 


where J C (P, k) = In k — 7 P 2 


is the current value objective function, and P fulfills the ODE initial value problem 


p(t j 2 

P{t) = k{t) - bPit) + 1 + y , Pi 0) = P 0 > 0. (lc) 

Here r, 7 , b > 0 are parameters, P = P(£) is the phosphorus contamination of the lake, which we want 
to keep low for ecological reasons, and k = k(t) is the phosphate load, for instance from fertilizers 
used by farmers, which farmers want high for economic reasons. The objective function consists of 
the concave increasing function lnfc, and the concave decreasing function — 7 P 2 ; b is the phosphorus 
degradation rate in the lake, and r is the discount rate. The discounted time integral in (la) is 
typical for economic (or socio-political) problems, where “profits now” weight more than mid or far 
future profits. More specifically, r often corresponds to a long-term investment rate. We focus on the 
parameter choice 

r — 0.03, 7 = 0.5, b G (0.5, 0.8) (primary bifurcation parameter), (2) 


and for the distributed case we shall additionally fix the diffusion parameter to D — 0.5. 

The max in (la) runs over all admissible controls k and (associated) states P; for k we can take 
the space C®([ 0, 00 ), R), and for P the space C£([0, 00 ), M). In fact, we naturally have k(t) > 0 for all 
t as J c (fc,P) —» —00 as k \ 0, and then (lc) implies that P(t) > 0 for all t > 0. On the other hand, 
see Remark 2.4 for comments on the case of state or control constraints. 

By Pontryagin’s Maximum Principle [PBGM62], see also, e.g., [GCF+08], an optimal solution 
(P*(-), fc*(-)) has to satisfy the first order optimality conditions 


k*(t) = argmaxP(P(t), fc, q(t), qo) for almost all t > 0, 
k 

with the local current value Hamiltonian function 

/ p2 

U(P, k, q, q 0 ) := q 0 J c {P, k) + q [ k - bP + ^—£>2 

The state P(-) and costate q(-) paths are solutions of the canonical system 1 

Pit) 2 

Pit) = d q H(P(t),k*(t),q(t)) = k*(t) - bP(t ) + 1 + y , 


with 


q(t) = rq(t ) - d P 'H{P{t),k*{t),q{t)) = 2y P(t) + q(t) ( r + b - 

P( 0) = Po > 0, 


2 Pit) 


(1 + Pit) 2 ) 2 >’ 


(3a) 

(3b) 


(4a) 

(4b) 


additionally satisfying the transversality condition 

lim e~ rt qit) = 0 if liminf P*(£) > 0. (5) 

t—>• oo t—yoo 


A solution (P(-), g(-)) of the canonical system (4) is called a canonical path , and a steady state of (4) 
is called a canonical steady state (CSS). Due to the strict concavity and continuous differentiability 
of the Hamiltonian function with respect to the control /c, and the absence of control constraints, the 
solution of (3a) is given by 


dkHiP{t), fc(t), q(t)) — 0 which yields k*(t) 


1 

qit)' 


( 6 ) 


1 It can be proved that the problem is normal, i.e. qo > 0, and hence w.l.o.g. qo = 1 can be assumed and is therefore 

subsequently omitted. 
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Consequently, for (P(-), </(•)) a canonical path, i.e., a solution of the canonical system, with a slight 
abuse of notation we also call (P, k ) with k = — 1/q a canonical path. In particular, if (P, g) is a CSS, 
so is (P, fc). Canonical paths yield candidates for optimal solutions, defined as follows: 

Definition 2.1. (P*(-), fc*(-, Po)) is called an optimal solution of (1) if for every admissible k(-) and 
associated P(-) we have 

./(ft,MO) < ./(ft,ft)) = V(Po). 

Then fc*(-,Po) is called an optimal control, P*(-) t/ie corresponding optimal (state) path, and 

m = k*(t, p 0 ) - bP(t ) + 1 ( 7 ) 

is called the optimal ODE. A constant solution (P*(-), fc*(-, Po)) = (P, k(P )) of (7) is called an optimal 
steady state (OSS). 

It turns out that the long-run behavior of an optimal solution (P*(-), &*(•)) can be characterized 
completely, see, e.g., [ Vag03]. Each optimal solution converges to an OSS, and depending on the 
parameters (4) can have I — 1, 2, 3 CSS (P, q)i, i = 1,...,/. 2 For simplicity omitting the non-generic 
case / = 2, if / = 1 then the unique CSS is a globally stable OSS, while for / = 3 two CSS are locally 
stable OSS, and the third is unstable. Here a OSS (P, q) is called globally (locally) stable if for each 
P(0) (in a neighborhood of P) the associated optimal path converges to (P,g); see [KW10, Kisll] for 
a detailed discussion. 

Setting u := (P, q) and letting u be a steady state of (4), the problem now is to compute a path, 
or “connecting orbit”, with P(0) = Po and lim u(t) — u. One standard approach, see, e.g. [LK80, 

t—► oo 

17, BPS01] and in particular [3CF+08, Chapter 7], is to treat (4) on a finite time interval 
[0,T] and to require u{T ) G W s (u ), where W s {u ) is the local stable manifold of u. In practice we 
approximate W s {u ) by the stable eigenspace E s (u ), and thus require 

u(T) G E s (u) and close to u. (8) 

To obtain a well defined two point boundary value problem we then need dim E s (u) = 1. 

More generally, if the state variable is an n-dimensional vector and thus the canonical system is a 
system of 2 n ODEs, for arbitrary Pq we need that u has the saddle point property, defined as follows. 

Definition 2.2 (Saddle Point Property). A CSS u G M 2n with 

n s := dim E s (u) — n (9) 

is called a CSS with the saddle point property (SPP). The number d (u):=n s — n is called the defect 
of u, and a CSS with defect u < 0 is called defective. 

For ODE problems like (4), given u with the SPP, and some initial state Po, canonical paths 
connecting Pq and u can now be computed using, e.g., OCMat. We now generalize this to distributed 
problems, and thus in §3 explain further details on that level. 

2.2 The shallow lake model with diffusion 

Following [ X08] we consider the shallow lake model with diffusion in a domain D C R d , d— 1,2, i.e., 

POO 

F(Po(0) :=maxJ(Po(.),fc(v)), J(P o (0, K, •)) := / e~ rt J ca (P(i), k(t)) dt, (10a) 

fc(v) Jo 

2 cf., e.g., the FCSS branches in Fig. 1 for our specific parameter choice. 
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where J ca (P(-, t), fc(-,t)) = f J c (P(x,t),k(x,t))dx (J C (P, fc) = lnfc — 7 P 2 as before) (10b) 

is the spatially averaged current value objective function, and P fulfills the initial boundary value 
problem 


dtP(x, t ) = k(x, t) - bP(x , t) + - + DAP(x, t ), (10c) 

i + P(x,ty 

d u P(x,t)\ gn = 0, P(ar,i)| t=0 = Po(^), a; G fl C R d , (lOd) 

where A = +... + d 2 rf , an d v the outer normal of O. We normalize J ca by the domain size \Q\ for 

easy comparison between the 0D, ID, and 2D cases, and, more generally, between different domains. 
We mostly focus on D = (—L,L) a real interval, but also give an outlook to 2D. In 2D the model is 
somewhat less intuitive, as a controlled phosphate dumping in the “middle” of the lake from farming 
appears difficult to motivate, and thus in 2D we rather think of (10) as a general pollution model. 
Instead of the periodic BC in ID in [ 1X08] we require Neumann (zero flux) boundary conditions (BC), 
which from a modeling point of view we find more natural. 

Introducing the costate q : D x (0, oo) —>> R^ and the (local current value) Hamiltonian 

p 2 

n = H(P, q, k) = J c (v , k) + q[k-bP + + DAP } , (11) 


by Pontryagin’s Maximum Principle for % = / 0 °° e rt P(t) d t with the spatial integral 

H(t) = j H(P(x,t),p(x,i),k(x,t)) dx, ( 12 ) 

Jn 

the canonical system for ( 10 ) becomes 

d t P(x,t ) = [d q H](x,t) = k(x,t ) - bP(x,t ) + 1 ppjppi + DAP(x,t), (13a) 

2 P(x, t ) 


d t q(x, t ) = rq(x, t ) - [dpH](x, t ) = 2-fP(x, t ) + q(x, t) [r + b - 


(1 + P(x,t) 2 ) 2 


d v P(x,t)\ dn = 0, d v q(x,t)\ m = 0, P(x,t)\ t=0 = P 0 (x), xett, 


— DAq(x, t), (13b) 
(13c) 


where k = argmax^ 'H(P, q, k ), which similar to ( 6 ) is obtained from 

d k U{P,q,k) = 0 & k(x,t ) = — - 1 — - (13d) 

q(x,t) 

The costate q also fulfills zero flux BC, and derivatives like dpT-L etc are taken variationally, i.e., for 
T~L. For instance, for <F(P, q) := qAP we have 4>(P, q) — f^qAPdx = f Q (Aq)P dx by GauB’ theorem, 
hence 5p$(P,q)[h\ = f (Aq)hdx, and by the Riesz representation theorem we identify 5p<h(P, q) and 
hence dp$(P, q) with the multiplier A q. 

Finally, we have the limiting intertemporal transversality condition (see Remark 2.4 below) 

lim e~ rt [ q(x, t)P(x, t)dx = 0. (14) 

Jn 

Analogous to Def. 2.1 we define 
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Definition 2.3. Let (P*(-, •), /c*(-, *, Pq)) be an optimal solution of problem (10), i.e. for every admis¬ 
sible k (•, •) and associated P(-, •) we have 


J(Po, *(■, *)) < W, •, P 0 )) = V(P 0 ). 


Then fc*(-,-,Po) is called a (distributed) optimal control, P*(-,-) is called the associated distributed 
optimal (state) path, and 


d t P(x,t) = fe*(x,t,Po(^)) - bP(x, t) + 


1 + P(x,t) 2 


+ DAP(x,t), d t/ P(x,t)\ dn = 0 (15) 


is called the optimal PDE. Again with a slight abuse of notation, (P*,k*) is also called an optimal 
path, and an optimal stationary solution (P(-), &;(•)) of (13) is called an OSS (optimal steady state). 
//P(-) = P then the optimal steady state is called a FOSS (flat optimal steady state), otherwise it is 
called a POSS (patterned optimal steady state). 

For a CSS u(-) we additionally introduce the acronyms FCSS for a flat canonical steady state, i.e 
u(-) = u, and PCSS (patterned canonical steady state) otherwise. Obviously, the FCSS correspond 
precisely to the OD CSS from § 2 . 1 . It was already indicated in [3X08] that (13) can additionally 
have PCSS arising from Turing like bifurcations. Thus, we first calculate bifurcation diagrams for 
(13), in ID and 2 D, recovering the up to 3 branches of FCSS from § 2 . 1 , and many branches of PCSS. 
Next, analogous to the OD case, we expect a solution of (10) to converge to some CSS. Thus, we only 
consider solutions u(-, •) of (13) with lim^oo ^(*, t) = &(•), where u(-) is a CSS. For u we then also need 
a version of the SPP. However, E s (u ) (and W s (u )) and E u (u ) (and W u (u )) are infinite dimensional. We 
circumvent this problem by requiring (8) and the saddle point property after a spatial discretization, 
which turns (13) into a (very large) systems of ODEs again. See App. A for further discussion. 

Remark 2.4. (a) A strict mathematical proof of Pontryagin’s Maximum Principle for diffusion pro¬ 
cesses over an infinite time horizon is still missing, specifically for the transversality condition (14), 
which also in [13X08] is discussed only on a heuristic base. Thus, at the moment we apply Pontrya- 
gin’s Maximum Principle in an ad hoc sense. We specifically assume, based on the results for the OD 
shallow lake model, that canonical paths converge to CSS, and therefore make no use of the “critical” 
transversality condition (14). In any case, a particular feature of the canonical system for diffusion 
processes is the anti-diffusion in the co-states, cf. (13b), which makes the canonical system ill-posed 
as an initial value problem. 

(b) For background on OC in a PDE setting see also for instance [ 4*610] and the references therein, 
or specifically [RZ99a, RZ99b, DHM 12 , ACKLT13] and [AAC11, Chapter5] for Pontryagin’s Maximum 
Principle for OC problems for semilinear parabolic state evolutions. However, these works are in a 
finite time horizon setting, and often the objective function is linear in the control and there are control 
constraints, e.g., k(x,t ) E K with some bounded interval K. Therefore k is not obtained from the 
analogue of (13d), but rather takes the values from dK , which is usually called bang-bang control. 
In, e.g., [Neu03] and [DL09], stationary spatial OC problems for a fishery model with (active) control 
constraints are considered, including numerical simulations, which correspond to our calculation of 
canonical steady states for our SLOC model. Here we do not (yet) consider explicit control or state 
constraints, and have an objective function strictly concave in the control, and thus we have a rather 
different setting than the above works. 

(c) We summarize that we do not aim at new theoretical results, but rather consider (13) after a 

spatial discretization as a (large) ODE problem, to which we apply the “connecting orbit method”. 
Importantly, this means that we take a broader perspective than aiming at computing just one optimal 
control, given an initial condition Pq, which without further information is an ill-posed problem anyway. 
Instead, our method aims to give a somewhat global picture by identifying the pertinent CSS and their 
respective domains of attraction. | 
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3 Numerical algorithm, and results 

The general idea is to use a method of lines discretization of (13), i.e., to approximate 

2 n 

u(x, t ) := (P(x, t),q(x, t)) = ^ x), (16) 

2=1 

where spans a subspace X n of the phase space X of (13), e.g., here X = [Tf 1 ^)] 2 , and 

(ui,u n +i), i = 1,... , n, are the expansion coefficients of P, </, respectively. This converts (13) into a 
(high dimensional) ODE 


u(t) =-G(u(t)), Ui(0) = u 0ii , i = 1 ,..., n, (17a) 

for the coefficient vector u = (ui)i= i v .., 2 n 5 where we have initial data for exactly half of the expansion 
coefficients. We choose a truncation time T and augment (17a) with the approximate transversality 
condition 

u(T) G E s (u), and || u(T) — u\\ small, (17b) 

where u is a steady state of (17a), and E s (u) is spanned by the eigenvectors of —d u G belonging to 
eigenvalues with negative real parts. As in OD we then need the SPP, i.e. 

dim E s (u) = dim E u {u), (18) 

to chose an arbitrary initial point in the state space. For a discussion and possible extension of the 
SPP to PDEs see App. A. 

Arguably, the simplest discretization for (16), at least in ID, is a finite difference (FD) scheme, 
which has the advantage that we can directly use OCMat for (17a),(17b), see [Gral5]. However, here 
we opt for a FEM ansatz for (16), using the setting of pde2path [UWR14, DRUW14] for two reasons: 

1. We want to consider (10) and hence (13) also on general 2D domains, and for more general 
models where the state variables may be vector valued functions already, see §4, again in ID or 
2D. In all these cases, FD and the coding of the respective spatially discrete systems may become 
rather inconvenient, while pde2path provides convenient interfaces precisely for such systems. 
Moreover, for more complicated systems adaptive meshes may become important, which are 
more easily handled in a FEM discretization, and are already an integral part of pde2path. 

2 . As explained above, the canonical system may have many stationary states; it is thus desirable 
to use a continuation and bifurcation package to conveniently find CSS. The goal then is to 
“seamlessly” link the setting of pde2path for stationary problems with BVP solvers for (17a). 

On the other hand, a drawback of spatial FEM discretizations is that the associated evolutionary 
problems have the natural form 

Mu{t) = -Ku(t ) + MF{u{t)) =: -G(u(t)), (19) 

where u corresponds to the nodal values, M, K G M 2nx2n are called the mass matrix and the stiffness 
matrix, respectively, and F : M? n -G M 2n is the nonlinearity. Thus, M and K are large but sparse; the 
signs in (19) comes from the convention that pde2path discretizes —A as K (positive definite). 
The occurrence of M on the left hand side of (19) means that it is not of the form (17a), and creates 
problems for the usage of standard BVP solvers. Of course, M is non-singular, and hence (19) can 
be rewritten as 

u = -M~ 1 Ku + F(u), (20) 

where for speed we can pre-calculate M _1 K. However, M -1 is no longer sparse, and already for 
intermediate n (n > 1000 , say) this results in slow computations and in particular memory problems 
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when using standard BVP solvers, which sort the Jacobian —M _1 K + d u F from each time-slice 
to,..., t m into a big Jacobian for the BVP problem. This can be alleviated by providing an approximate 
Jacobian J = — A + d u F, where A approximates M _1 K via lumping, i.e. we drop entries from M~ X K 
below a certain size S. Of course, there is a trade off between accuracy of J and the number of its 
nonzeros. 

We mainly experimented with the Matlab solvers bvp4c and similar, for instance the adaptions of 
bvp4c already implemented in OCMat, and the solver TOM [MS02, MSTC ] 3 . It turned out that TOM 
worked best in the “lumped” setting. Moreover, TOM was easy to modify in an ad hoc way to handle 
M on the left hand side, and a new official release of TOM is scheduled that also uses M [Mazl5]. 
For speedup it is advisable to avoid numerical differentiation and hence to pass a Jacobian function 
J=fjac(t,u) to TOM. This is generically very easy as pde2path in most cases provides a fast and 
easy way to assemble J. See [ Jecl5t] for implementation details. 

3.1 The algorithm for the computation of a path to a CSS satisfying the SPP 

To calculate canonical paths from a given state Po that connect to some CSS u with the SPP we want 
to solve the two-point BVP 


Mu[t) = -G(u(t)), t G (0, T) (21a) 

Pi{ 0 ) = P 0 ,i? i = 1 ,..., n, (n left BC), ( 21 b) 

\I !{u{T) — u) — 0 G (n right BC), (21c) 

where T G R nx2n encodes the projection onto the unstable eigenspace, i.e. ^(u — u) = 0 for u G E s (u), 
and where T is the chosen truncation time. The calculation of J/ at startup, which for large n turns 
out to be one of the bottlenecks of the algorithm, also gives a lower bound for the time scale T via 
T > _ R 1 e/ii , where /ii is the eigenvalue with largest negative real part, i.e., gives the slowest direction 

of the stable eigenspace of u. In our simulations we typically use T between 50 and 100. 

In general, a BVP solver needs a good initial guess of t i—^ u(t) to solve problem (21). Therefore 
we embed problem ( 21 ) into a family of problems replacing ( 21 b) by 

P(0) = aP 0 + (1 - a)P, a G [0,1], (21d) 

where we assume that for some a the solution is known: this holds for instance for a = 0 with the 
trivial solution u = u. We may then gradually increase <a, using the last solution as the new initial 
guess. This is implemented in the algorithm summarized in Table 1. 

There are more sophisticated variants of the simple continuation in Step 2 of Table 1 (some of 
which are implemented in OCMat), but the simple version in general works well for the problems we 
considered. Nevertheless, it may be that no solution of (21a), (21c) and (21d) is found for a > for 
some oio < 1, i.e., that the continuation to the intended initial point fails. In that case usually the 
BVP problem undergoes a fold bifurcation. We then use an adapted continuation step 2! that allows 
us to continue solutions around the fold. 

Remark 3.1. (a) For some applications it is useful to rescale the time t — Tr and hence consider 
Mu(t ) = TG(u(t ))) on the normalized time interval r G [0,1], which turns the truncation time T 
into a free parameter. This is for instance implemented in OCMat, but for simplicity not used here. 

(b) Similarly to the normalized normalized objective value J ca in (10b), in the bifurcation diagrams 
we use the normalized L 2 norm for comparison between different domains and space dimensions, i.e., 
henceforth, ||P ||2 := ||P||^ 2 /y / |i 2 |, and in the table in Fig. 1 we present averaged values, i.e., 

{p):= w\L p(x)dx ’ {k):= w\L k(x)ix - 

3 see also http://www.dm.uniba.it/~mazzia/bvp/index.html 


( 22 ) 




Step 0 (selection of u and implementation of ( 21 c)). Given u we solve the generalized adjoint 
EVP d u G(u) T & = AM$ for the eigenvalues A and (adjoint) eigenvectors 4>, which also 
gives the defect d(u). If d{u) = 0, then from $ E C 2nx2n we generate a real base of E u (u) 
which we sort into the matrix 4/ E R 77 ' x2r h 

Step 1 (selection of initial mesh and initial guess). To start the BVP solver we choose the initial 
guess u(t) = u on a suitable initial grid 0 = to < t\ < ... < t m = T. Typically, we choose 
m = 20 at startup, and afterwards TOM uses its own mesh-adaption strategy. 

Step 2 (solution and continuation). Using (21d) we try to increase a in small steps 6 to a = 1, 
in each step using the previous solution as the new initial guess, often trying 5 = 1/4. 
After thus having computed the first two solutions we may use a secant predictor for the 
subsequent steps. 

Step 2! (arclength continuation). If the continuation fails for a > ao with ay> < 1, then we use 
a pseudo-arclength continuation for a modified BVP, letting a be a free parameter. 


Table 1: The continuation-algorithm iscont (Initial State Continuation); Steps 0,1 are preparatory, 
Step 2 or 2! is repeated. See [Uecl5b] for implementation details, and for remarks on performance. 


To take the finite truncation time T into account we let 

—rT 

J(K),T) :=J(P 0 (-),k(-r)) + ^rJca(P(T),k(T)). (23) 

Obviously, for T - the last term can be made arbitrarily small, while for CSS it yields the exact 
r 

(discounted) objective value. In the following we drop the tilde in (23), and write, e.g., for the 
objective value of a CSS u, and, e.g., Jp^u for the CS-path which goes from P\ to u. | 


3.2 ID canonical steady states 

Recall that first we use pde2path to study the steady state problem for (17a), i.e., 

o =-w)- bPix)+ TM? +DAp(x) - (24a) 

0 = 2 cP{x) + q(x) (r + b- 7~~rP^] ~ DAq(x), (24b) 

V (1 + P(x) 2 ) J 

9vP(x)\ dQ = 0, d u q(x)\ dQ = 0. (24c) 

In ID we choose D = (—L, L ) x (— 5 y , S y ) with Neumann BC on all boundaries and small 5 y such that 
we can use just two grid-points in y — ^-direction and the solutions will be constant in y\ this we call 
a quasi ID setup. The steady states of the canonical system for the 0 D model (1) are FCSS of (24), 
and for easy reference we introduce the acronyms in Table 2 . The FCSS are of course independent of 
the domain, but to search for PCSS bifurcating from FCSS, the domain size 2 L should be close to a 
multiple of 2i r/fc c , where k c is the wave number of a Turing bifurcation. The parameters in (2), with 
b near 0.7, yield k c « 0.44 [BX08], and for simplicity we then choose L = 27r/0.44 « 14.28. 

Continuing the FSI branch in b we find a number of Turing like bifurcations to PCSS, and follow four 
of these; see Fig. la,b for the bifurcation diagram(s), and (c) for example solutions 4 . On the branches 
pl,p2, and p3 there are secondary bifurcations, not further considered here. For the subsequent 
examples we focus on b — 0.75 and b = 0.65, check the SPP (18) for all CSS and find it only to be 

4 The notation, e.g., pl/ptl6 follows the pde2path scheme, e.g.: continuation step 16 on the branch pi is stored in 
folder pi and file ptl6.mat. 
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name description 

FSM Flat State Muddy, the upper FCSS branch with a high phosphor load P 

FSI Flat State Intermediate, the upper half of the second FCSS branch, intermediate P 

FSC Flat State Clean, the lower half of the second FCSS branch, low P 

Table 2: Classification of the FCSS branches, see also Fig.l. The high P state is also called eutrophic, 
and our “muddy” refers to the fact that under eutrophic conditions there are a lot of algae and other 
organic matter in the lake, while under oligotrophic conditions (low P ) the water is much “cleaner”. 


fulfilled for the FSM, for the FSC, and for some points on the pi branch, e.g. at point 71 between the 
folds (see Fig. 1). 


(a) BD of CSS, (normalized) L 2 (b) BD, current values J c ^ c 

norm over b 


(c) example CSS 




p1/pt16 


p1/pt71 


1 

0.5 

0 





name 


(d) Characteristics of points in (a)-(c). 
(. P) ( k) J d name (P) 


(k) 


J 


d 


FSM/pt20 1.22 0.32 -63.11 0 

FSM/ptll 1.44 0.26 -79.28 0 

FSI/pt36 0.87 0.13 -79.47 -5 

FSC/ptl2 0.45 0.12 -72.95 0 


pl/ptl6 0.61 0.14 -74.83 -1 

pl/pt71 1.24 0.22 -78.93 0 

p2/ptl6 0.76 0.15 -76.70 -2 

p3/ptl9 1.02 0.17 -79.48 -3 


Figure 1: Basic bifurcation diagrams (a) and (b), example plots (c), and characteristic values of 
selected CSS (d). In (a), the blue and black lines represent the FCSS, and for instance the red line pi 
contains patterned CSS with one “interface” between high and low P. The numbered points on all 
these lines correspond to selected solutions plotted in (c), and characterized in (d). The small circles 
in (a) denote bifurcation points. The values J ca and J — J C a/r of the CSS are all negative, but this 
is merely a question of offset. 


3.3 ID canonical paths 

For b — 0.75, the only CSS is the FSM (FSM/pt20). It has the SPP, we can reach it from an arbitrary 
initial state Po> and thus it is a globally stable FOSS. Therefore, this regime is not very interesting, 
and we immediately turn to the case b = 0.65 with multiple CSS. 

For b — 0.65 seven CSS are marked in Fig. la, and characterized in the table, where only three 
satisfy the SPP. These are the FSC (which has the maximal value among these CSS), the FSM, and 
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the PCSS pl/pt71, subsequently denoted as ups. Next we numerically analyze which of these CSS 
belong to optimal paths. 


3.3.1 PCSS not satisfying the SPP 


From the analysis of non-distributed optimal control problems we know that steady states that do 
not satisfy the SPP can nevertheless be optimal. To illustrate that this is at least not typical in the 
SLOC model, we first compare the objective value of some (arbitrary chosen) CSS, e.g., the PCSS 
fips-(-) p3/ptl9, not satisfying the SPP, with that of the ^-dependent canonical paths rq(-, •) which 
connect to fq(-), i G {FSC, FSM, PS}. For the objective values we write J PS - for the CSS, and, e.g., 
Jp S -^ FSC for the canonical path which goes from P pg - to ufsc- The optimal solution for Pq = P P g- 
has to satisfy 


v(p PS ~) 


maxjjpg-, T 


p ps-^ fsc ’ Jp ps-^ F S M ’ J P PS -^ PS 


}• 


(25) 


The canonical paths are given in Fig. 2a-c, while (d) presents some norms along the path in (a), which 
show that and how fast u(t) (including the co-states q) converges to u. In all cases we find without 
problems canonical paths to both FCSS and the PCSS; in particular the path to FSM is rather quick. 
For the objective values we find, up to 2 significant digits, 


Tps— 79.48 < Tpg-^ FSC — 78.24 < T PS -^ PS — 78.19 < T PS -^ FSM — 77.5. (26) 

Thus, the optimal path is the one converging to ufsm- Repeating these steps for every PCSS not 
satisfying the SPP we find that these are always dominated by paths converging to one of the FCSS. 
Therefore, only ufsc, ^fsm and tips remain as candidates for OSS. 

Before we turn to determining OSS, we briefly discuss the canonical paths in Fig. 2. We focus on 
(a), but similar remarks apply to (b,c) as well. In (a), the initial P(-,0) is above the target Pfsc? so 
naively we may expect that the control k should start below the target fcpsc an d slowly increase to 
fcpsc- However, such a control would not be optimal. Instead, k is initially similar to /c PS -, and in 
particular fc(-, 0) is large where Po(-) is already large. Only after a short transient k drops below fepsc 
and then behaves as expected. 

At first sight, this startup behavior of k may appear rather counter-intuitive. However, the reason 
is that we do not want to drive the system to ufsc “as quickly as possible”, which essentially would 
amount to choosing k “as small as possible” at startup. Instead, we want to maximize J, and for this 
it pays off to have, for a short transient, k large near the maxima of P(-,0). 

To illustrate this point, in Fig. 2(e) we choose the naive control k(t ) = fepsc for a ll t and numerically 
integrate the initial value problem (10c). While this does take us to ufsc, the first observation is that 
this needs a rather long time. Secondly, for the value of this solution we obtain J = —80.1, which is 
even worse than the starting CSS with J PS - = — 79.48. 


3.3.2 Determining optimal steady states 

We return to the question whether one or more of the CSS at b — 0.65 with the SPP from §3.3.1 
are optimal, and proceed in three steps. First we search for a canonical path starting at P P s and 
connecting to ftpsc- In the second step we repeat that for ufsm, and in the last step we check if both 
or only one of the FCSS are optimal. 5 

Paths between Pps and ufsc — a Skiba candidate. Using iscont to get a canonical path 
starting at P P g and converging to the FSC it turns out that the continuation (21d), i.e., 

Pa( 0) := aPps + (1 - a)PpsCi (27) 

5 The second step reveals that the first step is superfluous, but this we do not know a priori. 
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(a) canonical path from p3/pt!9 to FSC 


(b) canonical path from p3/pt!9 to FSM 



Figure 2: Canonical paths from the (state values of) PCSS - p3/ptl9 to FSC (a), FSM (b), and PS 
(c) at b — 0.65, and typical path diagnostics (d). For comparison, (e) shows the solution P of the 
initial value problem (10c) with P(0) as in (a) and the externally chosen control k(t) = &fsc for all t. 


yields a fold around a « 0.6, see the blue curve in Fig. 3a), and that no canonical path connecting 
Pps to upsc exists. Instead, multiple solutions that converge to upsc exist for initial distributions of 
the form (21d) with a E [0.6, 0.71]; two examples are shown in Fig. 3b, and their diagnostics in (c). 

Similarly, trying to continue to a path that connects Ppsc and tips yields a fold (green curve in 
(a)), and no such path exists. However, the solutions returned during the continuation process allow 
us to determine and compare the respective objective values. This yields that there exists a specific 
initial distribution where the objective values are equal, given by the intersection of the green and blue 
curves in Fig. 3a. Thus, from an economic point of view both solutions are equal. This suggests that 
Ps is a Skiba or indifference threshold point (distribution) 6 , well known from non-distributed optimal 
control problems, see, e.g., [ Vag03, KW10]. The paths u starting at Ps (red curve) are depicted in 
Fig. 3d: P(-,0) is the same for both solutions, but the controls fc(*,0) are different. In any case, to 
assure that these solutions are optimal we have to prove that no other dominating solution exists. 
Thus in a last step we calculate the objective values of the paths converging to uysm- 

From P PS to ^fsm* Here the continuation is successful and we find a path connecting Pps to ufsm- 
Comparing the objective values reveals that the PCSS is dominated by the solution converging to the 
FSM, see Fig. 3e and 3f. Thus, the PCSS is ruled out as an optimal steady state, and therefore we a 
posteriori identify Ps as only a Skiba candidate as it does not separate two optimal steady states. 

A Skiba manifold between PpSM and Ppsc* It is well known that in 0D the FSC and FSM are 
only locally stable with regions of attractions separated by a Skiba manifold (parametrized by, e.g., 
b ) of homogeneous solutions, [KW10], see Fig. 3g for our case b = 0.65. Of course, this also yields a 

6 however, taking into account also the FSM solution, below we identify Ps as only a Skiba candidate 
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(c) Diagn. for left and right paths in (b) 



(d) A Skiba candidate between Pps and FSC (e) p, a for FSM target 


(f) Dominating path to FSM 



(h) Homog. Skiba point and 
(g) OD Shallow Lake paths to FSC and FSM 



(i) A patterned Skiba point between FSC and FSM 



Figure 3: Canonical paths to various CSS for b — 0.65, and illustration of some Skiba points; see text 
for details. 


homogeneous Skiba distribution in ID, see Fig. 3h. More generally, we may expect the domains of 
attraction of the FSC and the FSM to be separated by an (in the continuum limit infinite dimensional) 
Skiba manifold Mg, for fixed b. 

A continuation process, analogous to the non-distributed case, see [Gral2], could be used to ap¬ 
proximate this manifold Mg. However, to find a non homogeneous example point on that manifold, 
here we can readily combine Fig. 3a and 3e, to find a Skiba distribution of the form 

Pskiba = o'Pfsc + (1 - <a)Pps; (28) 

see Fig. 3i for paths to the FSC and the FSM yielding the same J — —76.3. 

Summary for ID. The picture that emerges is as follows: for b > bf Q \d ~ 0.727 the FSM as the only 
CSS is the globally stable FOSS, while for b < b^ there exist multiple CSS. Specifically for b — 0.65, 
fiFSC 5 f^FSM and one of the PCSS have the saddle point property, but only ufsc and ^fsm are optimal, 
and in particular no POSS exists. The FOSS FSC and FSM are separated by a (presumably rather 
complicated) Skiba manifold Mg, and Fig. 3h and 3i show just two examples of points on Mg. 
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3.4 Outlook: 2D results 


As a 2D example we consider (13) on the domain D = (—A, L ) x (—L/2, L/2 ), L = 27t/ 0.44 as before, 
with a rather coarse mesh of 40 x 20 points, hence approximately 1600 DoF. The FCSS branches are 
of course the same as in ID (or 0D), and again at the end of the FSI branch we find a number of 
Turing like bifurcations. In Fig. 4(a), (b) we only present the “new” patterned branches, i.e., those 
with a genuine x and y dependence. 


(a) Two 2D PCSS branches 


(b) some example plots 



P at h2/pt8 


Pat h2/pt17 




0.8 0 

0.6 -5 



-10 


10 


-10 


10 


P at h3/pt20 


I:, IL 10: 


-10 


10 


k at h2/pt8 k at h2/pt17 



k at h3/pt20 



-10 0 10 


Figure 4: New patterned branches in 2D; (x,y) G (—L,L) x (—f, §). 


These new bifurcating PCSS again do not fulfill the SPP. As an example for a canonical path, 
in Fig. 5 we present snapshots from a path from P of the PCSS h2/ptl7 to the FSC (see also 
www.staff.uni-oldenburg.de/hannes.uecker/pde2path for the movie), which yields a higher J 
than the PCSS, i.e., 

J(PCSS) = -77.53 < J(PCSS -> FSC) = -76.23 < J(FSC) = -72.97. (29) 

Thus, this PCSS is not optimal, and neither is any other one we checked. Using the methods from 
§3.2 it is now of course also possible to find points with a genuine x and y dependence on the Skiba 
manifold separating FSC and FSM, but here we skip this presentation. 



Figure 5: Solutions on the canonical path to FSC 


The behavior and economic interpretation of the path from the PCSS to the FSC in Fig. 5 is 
rather similar to the convergence to the FSC in Fig. 2a. After a short transient the optimal strategy 
is to give a high phosphate load k where P is below the limit value Ppsc (south-west and north-east 
corners of the domain), but initially there also is a high k at high P values (north-west and south-east 
corners). 
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4 Discussion 


We have presented a numerical framework to treat infinite time horizon spatially distributed optimal 
control problems. First we derive the canonical PDE systems, which we then discretize in space and 
thus approximate by large systems of ODEs. For these we can resort to the theory and experience with 
non-distributed optimal control problems. Thus our results are intrinsically numerical; however, we 
believe that they can help to develop the theoretical concepts for distributed optimal control problems, 
following Oliver Heavisides’ saying “Mathematics is an experimental science, and definitions do not 
come first, but later on. ” 

From the economic point of view, the computation of canonical paths to the FOSS yields nontrivial 
and interesting results. Even more interesting would be locally stable POSS, but there is strong 
evidence that these do not exist for the shallow lake model (10), at least in the parameter regimes we 
considered so far 7 . On the other hand, in [Uecl5a] we use our method to study the vegetation system 
from [ 1X10], and find that POSS dominate in large parameter regimes. We believe that the same 
happens in many other important systems, and a number of further investigations in this direction are 
under way. Natural candidates, i.e., systems with natural objective functions and controls, are related 
vegetation systems as in [SZvHMOl, ZKY+13], fishery models as in [Neu03, Gral2], or “crimo-taxis” 
systems as in [SBB10]. 


A SPP for PDEs 


In this appendix we discuss the SPP (Def. 2 . 2 ) in a somewhat more general situation, tailored to 
canonical systems coming from spatial discretizations of PDEs. Let u — (p, q) G M? N be a stationary 
state of a (non-distributed) canonical system of the form 


s(|i=^= 


f(p, <?) 
rq - H p (p,q))J ’ 


(30) 


where / = H q , and let J = D u F{u) be the Jacobian at u. In [GCF + 08, Thm 7.10] it is explained that 
the eigenvalues of J are symmetric around r/ 2, i.e., that there exist N complex numbers ^ such that 


<7(J) = {£±6:* = 1,...,JV}. 
In detail, since det(J — £) = det(J r — (£ — §)) where 


j j_L- H m- r 2 H « 

J r • 0 - tt _tt i r |) 

Z \ 11 pp 1J -pq \ 2; 


L PQ 


(31) 


(32) 


we have that | ^ G C is an eigenvalue of J if and only if ^ is an eigenvalue of J r . But J r has the 

f A B \ 

with symmetric matrices B,Cg M iVxAr , and as a consequence the eigenvalues of 


structure 


C -A 


J r are & = i 

Now consider the distributed canonical system 


= 1,...,W 


d t 


p(x, t ) 


= F(p(x,t),q(x,t)) + 


1 \q(x,t) 

where D € W NxN is a diffusion matrix, i.e., positive definite. Let 


( DAp(x,t) 
\-DAq(x,t) 


—u(t) = G(u(t)), u € 


b2nN 


(33) 


(34) 


7 see however “Scenario 2” in [Gral5] for some locally stable POSS 
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be the associated spatially discretized system with n spatial points, where u = (p(aq), ... ,p(x n ), 
q(x i),... ,q(x n )) G R 2n7V , and let u G M? nN be a steady state of (34). Then J = D u G(u ) has the 
structure J — —K + Ji oca i, where Ji OC ai has the block structure 


Aocal — 


( 

0 


_ TTl 
11 pp 


0 


tt2 

pq 


_ tt2 

12 pp 


0 

0 




H: 


m 


0 

0 

-H, 


TTl 

IP 




0 

0 

r - H 2 
1 qp 


pp 


0 

0 


h: 


qq 

0 

0 


r - H n 

' qp/ 


composed of local matrices H pq := H pq (xj) H pq (p(xj ), q(xj )), H qq , H pp G R NxN , and K — 


(35) 


L 0 


with L € R nN coming from the discretization of DA. The notation K of course reflects the FEM 
background of the present paper, but the same structure ^ ) occurs for any discretization, in 


x 0 -L, 

any space dimension, and for any D not necessarily diagonal, i.e., containing cross diffusion. 

It follows that again | ^ is an eigenvalue of J if and only if £ is an eigenvalue of J r := J — 

where J r has the structure 


j r — s y mmetric 5,c g R nN . 

Applying [GCF+08, Lemma B. 2 , Lemma B.3] we obtain 

Theorem A.l. Let u be a steady state of the spatially discretized distributed system (34), and let J 
be the associated Jacobian. Then there exist ^ G C, i = 1,..., nN, such that 

a(J) = { r -±^:i = l,...,nNy (36) 

As a consequence, dim E s (u) < Nn, and the only candidates u for right BC in ( 21 c) are those with 
the SPP. As a corollary we find a property that, on the discretized level, is equivalent to the SPP. 

Corollary A.2. Let u € R 2nN be an equilibrium of the spatially discretized distributed system (34) 
and r > 0. Then u satisfies the SPP iff every eigenvalue £ of the according Jacobian J(u) satisfies 

H Re€ - P > 5* (37) 

Remark A.3. Theorem A.l is formulated on the discretized level, and one might ask how it ultimately 
relates to the PDE. As a first step one can ask: Let a steady state u G M? nN of (34) be an approximation 
of a PDE steady state (p,q) G X for (33), with X C {(p, q) : D R 2iV } some function space, e.g., 
X — [iL 1 (D)] 2Ar . If u G M? nN is an approximation of (p,q) on a finer mesh h > n, or just a different 
mesh, do we have 


n N — dim E s (u) = n N — dim E s (u) ? (38) 

We do not want to go into the details here, but if E c (u) = 0, i.e., cr(J) H {Re£ = 0} = 0, then 
(38) is true, for large enough n, h. Given some u, this can be easily tested numerically, and it is also 
clear from an analytical point of view. Refining u to u we essentially add high frequency modes to 
the FEM (or finite difference) mesh. These introduce the same number of additional eigenvalues at 
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large positive and negative £ for the linearization J, because Jf(p,q) : X —>► X is relatively compact 
with respect to the Laplacian, i.e., w.r.t. (p,q) ^ (DAp,—DAq). On the other hand, the small 
eigenvalues pi, |/i^| < i? for some fixed R, are only slightly perturbed, i.e., | fi{ — Pi\ < C\\u* — u ||, 
where u* is suitably defined, for instance by interpolating u to the mesh of u. But || u* — ?i|| —> 0 as 
n,h —>> oo, which yields (38). 

To make this rigorous, we need to define appropriate function spaces and study the approximation 
properties of the spatial discretization. This is easy, as the stationary problem for (33) can be written 
as an elliptic system, and hence (p, q ) is arbitrary smooth, but we omit the details here. 

In fact, (37) can also be formulated on the PDE level and might therefore replace the SPP from 
Def. 2.2 for spatially distributed models. However, we also postpone an in depth analysis of this to 
future work. 
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